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Abstract 

We study a version of the Tangled Nature model of evolutionary ecology rede- 
fined in a phenotype space where mutants have properties correlated to their par- 
ents. The model has individual-based dynamics whilst incorporating species scale 
competitive constraints and a system scale resource constraint. Multiple species 
arise that coexist in a species interaction network with evolving global properties. 
Both the mean interaction strength and the network connectance increase rela- 
tive to the null system as mutualism becomes more extensive. From a study of 
the dependence of average degree on the resource level we extract the diversity- 
connectance relationship which conforms to the hyperbolic form seen in field data. 
This is adjudged to arise as a consequence of the evolutionary pressure to achieve 
positive interactions. The network degree distributions conform more strongly to 
exponential than to the null binomial distributions in all cases. This effect is be- 
lieved to be caused by correlations in the reproductive process. We also study 
how resource availability influences the phenotypical lifetime distribution which is 
approximately of power law form. We observe that the mean lifetime is inversely 
related to the resource level. 

1 Introduction 



A functioning ecosystem relies on various environmental resources to sustain its com- 
ponent species with some of these resources being universally required for survival. The 
limitation of supply rate or overall quantity of such a resource, acts as a restriction to 
the development of an ecosystem and systematic properties are expected to depend on 
the level at which this occurs (Loeuille & Loreau, 2003). It is known from field data 
that the diversity of species changes across energy and water gradients such as is found 
in latitudinal variations (Hawkins et al., 2003) (Bonn et al., 2004). In general, heteroge- 
neous environments will support greater species numbers at higher resource availabilities 
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with the relationship taking a monotonically increasing form (Currie, 1991) (Waide et 
al., 1999) (Bonn et al., 2004). This difference in diversity and the knock-on effect of 
elevated competition could conceivably lead to a change in the topology and dynamics of 
the species interaction network. These networks and their associated interaction matri- 
ces have been studied and are considered to be significant in adjudging the stability and 
permanence of ecosystems (May, 1974) (Tregonning & Roberts, 1979). 

There has been a great deal of research investigating the short term effects of resource 
variation on ecosystems (Tilman, 1982) (Hulot et al., 2000), but very little on the im- 
pact of resource levels over evolutionary timescales. Also, recent work has focused on 
the dynamical coevolution of species populations (Tokita & Yasutomi, 2003) (Drossel et 
al., 2004) (McKane, 2004) (Coppex et al., 2004) but these approaches use continuum de- 
scriptions whilst evolution occurs via dynamics associated with discrete individuals. To 
address these issues accordingly we present a model of evolutionary ecology that functions 
at the level of the individual whilst incorporating species scale competitive constraints 
and a system scale resource constraint. 

Our work follows from the Tangled Nature model of evolutionary ecology (Christensen 
et al., 2002) (Hall et al., 2002) (Collobiano et al., 2003) (Anderson & Jensen, 2005), which 
demonstrated species formation as a result of individual-based dynamics. These species 
were then studied in terms of their interaction networks and other more dynamical fea- 
tures such as extinctions, lifetimes and stabilities. Recently, the model was simplified 
by restricting the reproductive process to non-overlapping generations (Rikvold & Zia, 
2003) (Zia & Rikvold, 2004). This allowed both longer timescales to be simulated and a 
deeper analytical treatment to be made available. Modifications made to the description 
of individuals have now allowed us to include a gradual species evolution within a signif- 
icantly larger phenotype space. Correlations introduced into the structure of this space 
mean that inheritance from parent to mutated offspring is available and also mathemat- 
ically quantifiable which was not the case in the previous model. Attempts have been 
made to impose correlations onto the hypercubic genotype structure of the earlier model 
(Sevim & Rikvold, 2005) but with only limited effect. Here we use a different procedure 
to construct a correlated phenotype space which allows us to consider much larger spaces. 
We find that correlations in these huge phenotype spaces have a significant effect. 

Our paper is arranged as follows. In the first part of section 2, we describe how 
the system members are represented in terms of a hypothetical phenotype space. This 
representation is used to provide a fully determined, correlated interaction-space. As the 
processes involved in this are somewhat laborious we have consigned the wider details to 
the appendix, whilst giving an overview in the main text. In the second part of section 
2, we describe the update rules for the system and how the reproductive probabilities are 
determined from the species interaction network. We approach the dynamical evolution 
in a stochastic discrete manner as we see it as appropriate for a discrete-entity system 
such as this. Statistical fluctuations are important as population numbers may become 
low and coupled dynamics across an interaction network could be heavily dependent upon 
them. 

Results regarding the issues of diversity, lifetime distributions and network properties 
are presented in section 3, and a discussion of these along with future directions are set 
out in section 4. 
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2 Model 



2.1 Interactions 

Individuals are represented by vectors, T a = (Tf, T^,..., T"), in a phenotype space 
of dimension L = 16, with specific individuals denoted by Greek lettering a, (3, ... = 
1,2, ...N(t), where N(t) is the population at time, t. Each trait Tf may take an integer 
coordinate in the range [0,99999] that is periodically bounded allowing the points Tf = 
and T" = 99999 to be contiguous. The phenotype coordinates themselves are arbitrary 
and are not intended to represent any kind of quantitative scale so the periodic boundary 
is introduced to maintain their arbitrariness. A large integer range is used to emulate the 
continuous nature of many real organism traits and as a result the extant phenotypes are 
capable of gradual adaption during the evolutionary process. 

The co-evolution of the extant phenotypes is primarily controlled by the interactions 
they have with one another. These are represented by J alS = J(T a , T^), the interaction 
strength effected upon a by [3, which is independent of J(T /3 ,T Q? ). As the phenotype 
space is a closed set of 100000 16 = 10 80 possible states the interaction space is a closed 
set of 10 80 xl0 80 = 10 160 possible pairwise interactions that exist in potentia. This set of 
all possible phenotype interactions is constructed with the property that only a subset of 
them are non-zero quantities. This proportion or connectance, 6q is defined at the outset 
and represents the fact that a specified phenotype will interact only with a subset of 
all other conceivable phenotypes. The resulting set of non-zero interactions are assigned 
strength values that are normally distributed. 

A system created with a random group of phenotypes will, on average, conform to the 
above properties (ie. connectance, 6 = 6q, plus interaction strengths that are normally 
distributed). But, once the evolutionary process is underway the phenotypes existing at 
later times may have interaction properties that deviate from those of the random null 
system. The fact that we pre-define the properties of the closed interaction space means 
that the evolving system properties can be compared quantitatively to the null system. 

The phenotype space is constructed in a manner that gives a quantifiable correlation, 
C a/3 between separate points. This correlation measure decays approximately exponen- 
tially with separation in the phenotype space and has an effect on the inheritance of both 
the interaction set distribution (the subset of non-zero interactions) and the strength 
measures themselves. The consequence is that any mutated offspring will have a similar 
interaction set and strengths to those of the parent, with the level of similarity being 
dependent upon the mutation length. 

The requirement of a large, correlated, pre-defined interaction space means using a 
deterministic procedure for producing the interaction strengths. Any interaction between 
defined phenotypes must be time- invariant so as there are many more possible interactions 
than can be stored computationally we used a deterministic method for acquiring the 
values from any given pair of phenotype vectors. The exact method, whilst not difficult, 
is somewhat convoluted so is explained separately in the appendix along with the details 
of how the correlations are imposed. 

Intra-specific competition is deemed a necessary part of the dynamics to allow diversity 
and is incorporated by setting all self-interactions to a negative value that is constant and 
independent of the phenotype vector. This property extends via the imposed correlations 
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to phenotypes of similar but distinct vectors thus allowing for the localised but distributed 
nature of a species description. 



2.2 Dynamical update 

We initiate the system with a set of individuals which are assigned random phenotype 
vectors. At each subsequent timestep an individual is selected for annihilation with 
probability Pkm = 0.2 whereby it is removed from the system and the single resource 
unit associated with it is returned to the resource bath, R(t). During the same timestep, 
another individual is randomly selected to reproduce with probability P rep ro- This value 
is determined via use of the weight function, 

T N(t) r p N(t) 

ff(T " t)=0i ^^- 02 £ c °"- a w (1) 

This sum is monotonically mapped to the interval [0,1], appropriate for a probability 
measure, through use of the function, 



exp[H(T a , t)] 
l + exp[H(T a ,t)]' 

The first term in Eq.(l), represents the summed effect of interactions with other 
system members. Although it is made over all individuals at time, t, many of the con- 
tributions will be zero as the two phenotypes in question may not interact. This is an 
inherent feature of the complete interaction set which permits only a prescribed propor- 
tion, 9 , of all possible phenotype couplings to exist (as evolution occurs though, the 
connectance of the extant system, 6 may deviate from, 6> ). 

The denominator represents the number of individuals that can be described as be- 
longing to the same species group as T a . This damping effect is used to prevent diver- 
gence, but can be seen as realistic if we consider the sharing of the interaction effects 
amongst members of the same species. An example would be the predatory effect of a 
single lion amongst a herd of wildebeest. As the lion feeds up to a maximum rate the 
effect of the lion on an individual in a small herd is more pronounced than the effect in 
a large herd. Essentially, the probability of a specific lion consuming a specific wilde- 
beest decreases in relation to the size of the herd. This damping effect is by no means 
ubiquitous but it is nevertheless a feature of many interactions. 

Phenotypes that are in close proximity in the phenotype space are distinct but cor- 
related in a manner that can be described as a species set, so this must be accounted for 
here. This is achieved by using C a/3 G [0, 1], which is a measure of correlation between a 
and p. A more detailed explanation is given in the appendix. 

The second term of Eq.(l), represents conspecific competition so utilises the same 
correlation measure as the denominator in the first term. This negative contribution is 
essential to allow diversity in the system and generalises the multitude of factors other 
than the conserved resource that members of the same species compete over. 

The third term of Eq.(l), represents the general system competition for the conserved 
resource. It represents the number of individuals, N(t), that are competing per unit of 
available resource, R(t). As the condition for conservation, R(t) + N(t) = constant must 
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hold, the numerator and denominator are dependent. The resource range investigated 
here spans from 1000 to 30000 units. 

The parameters, a x = 0.5, a2 = 0.01, 03 = 0.2 are the selection, conspecific compe- 
tition and resource competition parameters respectively. These values were subjectively 
selected to allow significant diversities to be achieved whilst incurring the dynamics to 
be strongly dependent upon inter-specific interactions. 

Upon successful reproduction, the individual produces a single offspring that assumes 
one unit of resource from the resource bath. As the bath depletes the probability of 
reproduction for any phenotype is reduced due to the third term in Eq.(l). The offspring 
phenotype is identical to that of the parent (clonal) unless, with probability P mu t = 
0.0002, a mutation occurs. This is effected by shifting a randomly selected trait value by 
a random amount ip, that is gaussian distributed with /i — and a = £, where £ = 250 
is the correlation length of the phenotype space. Any mutated offspring will have an 
interaction set that is similar but not identical to its parent in accordance with the 
correlation. We have elected to use a low mutation rate with a relatively high mutation 
distance to allow the evolution to occur efficiently enough for tractable simulation times. 

3 Results 
3.1 Diversity 

By initialising the system with 9 = we can generate a neutral evolution in which 
the population increases whilst diffusing uniformly through the phenotype space. The 
intra-specific competition term causes this diffusion by forcing the system members to 
be as little correlated as possible. As a result, the diversity grows to large values whilst 
localised phenotype populations remain low. The incorporation of a non-zero #0 breaks 
this symmetry allowing phenotypes to counteract the competitive constraint with positive 
interactions and so accumulate localised populations. These phenotypes are distributed 
as highly populated single sites surrounded by a sparse cloud of mutants that derive 
primarily from the central 'wild type'. 

Fig.(l), shows a section of the time evolution of the extant species in a single run 
of one million generations. The visual representation of this is as a projection of the 
populated points of the 16 dimensional phenotype space onto a single trait. It is clear to 
see that the evolutionary process creates a system that is far from diffuse with a small 
set of phenotypes interacting in a manner that precludes easy invasion by mutants. Of 
course, there are successful invasions that amount to gradual evolution of a species, or 
even speciations, but the relative permanance of species is seen as significant. This is 
because a new mutant phenotype will have an advantage over the parent due to the 
relative weakness of its intra-specific competition term, so a continual invasion of species 
could easily be expected. The phenotype distribution localises at points rather than 
following a diffusive process and does so to quite an extreme. There is nothing to prevent 
the diversity from expanding with species achieving smaller populations but this state 
might struggle to persist. It is likely that the stochasticity of the dynamics would incur 
a greater extinction rate for species of such sparse numbers thus reducing the diversity. 
This is one reason proposed to explain why productivity-diversity relationships having 
increasing functional forms at low productivity ranges (Preston, 1962) (Abrams, 1995). 
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Figure 1: An occupation plot of a single run for a system with R= 10000. For each timeslice 
a point appears where a phenotype is in existence but as the full space is in 16 dimensions a 
projection onto a single trait is used. 



The diversity varies considerably both in time and across actualisations. This is most 
apparent for higher resource levels where the standard deviations of the diversity become 
comparable to the means. Regardless of this spread, for the range investigated here the 
mean species diversity increases with respect to total resource availability in a monotonic 
fashion, Fig. (2). This relationship has been produced in a species level trophic network 
model (McKane, 2004) and is empirically found in large scale systems with heterogeneous 
environments (Currie, 1991) (Waide et ah, 1999) (Bonn et ah, 2004). Although unimodal 
relationships are expected for localised ecosystems where diversity is more dependent 
upon fewer limiting factors this model incorporates aspects that are best compared to 
heterogeneous larger scale systems. At this ecological scale a monotonically increasing 
diversity can in some cases be ascribed to the effects of intra-specific density dependence 
(Abrams, 1983) (Vance, 1984) (Abrams, 1995). Resource increases allow species to grow 
in population but other factors more unique to that species' niche restrict this growth prior 
to the resource depletion becoming a limiting factor. The consequence is that resource is 
more freely available for species holding dissimilar niches that would be excluded at low 
resource due to their inferior ability to procure it. This model represents such systems 
as the intra-specific competition is the dominant restrictive term in Eq.(l), for a high 
species population. 

3.2 Lifetimes and extinction rates 

Statistical analyses of fossil record data have often alluded to power law forms, P(s) ~ s~ a . 
in the distributions of the various quantities involved with species extinctions. When s 
represents species lifetimes or extinction event sizes an exponent of a ~ 2 has been sug- 
gested but the analyses have been criticised and power law forms are not readily accepted 
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Figure 2: Plot of mean species diversity in relation to resource availability. Error bars represent 
the standard error. 



(Newman & Palmer, 1999) (Drossel, 2001). The lifetime distributions produced in this 
model are clearly not of this form although they do loosely follow a power law with a 
comparable exponent Fig. (3). Other models have displayed similar off-power law or non- 
power law behaviour (Standish, 1996) (Christensen et al., 2002) (Chowdhury et al., 2003) 
(Chowdhury & Stauffer, 2003a) (and also (Chowdhury & Stauffer, 2003b) in reference to 
(Sole & Manrubia, 1996)) giving credibility to the notion that these distributions may not 
actually be scale free in the real system. The power law form seen in real macroecological 
data is tenuous as a result of the difficulty in appropriating quality data sets, but it may 
be the case that the actual distributions are really deviations from this form. It is hard 
to draw conclusions from models such as these as they are obviously simplifications of 
a highly complex, spatially-extended, multiscaled system but there is no a priori reason 
why power law behaviour should be expected anyhow. 

Although the structural form of the distributions for different resource levels are sim- 
ilar there is a gradual decrease in the level of species permanence as resource increases. 
This effect can be seen in Fig. (3), as an increase in mid logarithmic-range lifetime prob- 
ability accompanied by a sharper declining tail. It is a minimal but consistent trend, as 
can be observed in the relationship between mean lifetime and resource level in Fig. (4). 
We can attribute this behaviour to the difference in mutant offspring production rates 
incurred at different resource levels. As the mean population increases with increasing 
resource availability there is a higher magnitude of successful reproduction events in each 
generational time step. So, given a constant mutation probability, the number of mu- 
tated offspring produced in this period increases in proportion to the population. The 
permanence of individual species and the stability of the interaction network are affected 
by this elevated mutant production rate leading to a greater rate of extinction and speci- 
ation. If we plot the mean lifetime of a given resource level against the mean population 
of that level Fig. (5) a linear dependence can be seen that supports the explanation above. 
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Figure 3: Plot of the lifetime distributions over all resource availability levels. In the approxi- 
mate range s = [W 2 , 10 3 ] the probability of a given lifetime increases with resource availability, 
whilst the reverse is true for the approximate range s = [10 4 , 10 5 ]. This indicates longer lived 
species at lower resources. The dotted line represents the often proposed power law relationship 
with exponent, a = 2. 
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Figure 4: Plot of mean lifetime against resource availability. 
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Figure 5: Plot of mean lifetime against mean population for each resource level. 



To justify this further we have run simulations with variable mutation probabilities that 
force the absolute mutant production rate to remain constant. The lifetime distributions 
produced from these runs show no dependence on resource availability which is consistent 
with our reasoning. 



3.3 Network properties 

An important feature of a coevolving ecosystem is the species interaction network, par- 
ticularly as it is thought to have implications for the stability of the ecosystem and the 
permanence of its members. Low diameter networks with scale invariant degree distri- 
butions have often been expected as many abiotic networks display such properties and 
there are theoretical reasons why they may augment stability (Albert & Barabasi, 2002) 
(Dorogovtsev & Mendes, 2002). In reality, field data has shown the existence of a variety 
of distributions including power law, truncated power law and exponential (Dunne et al., 
2002) (Jordano et al., 2003), but these still differ from the binomial distributions found 
in finite random networks. 

The degree distributions produced by our model take an exponential form which 
is similar, at low degrees, to the binomial distribution of the same connectance. The 
difference becomes apparent at higher degrees where the binomial decays much more 
rapidly than the exponential. 

Fig. (6), shows the degree distribution for the realisations with diversity of 29 (R = 
30000) along with a binomial distributed plot using the same connectance. It is clear 
that the evolved system has a distribution with a much longer tail than the binomial. 

The deviation of the degree distribution away from the binomial can be attributed 
to the inheritance-based development of the network structure. Simulating the evolution 
using uncorrelated mutation leads to binomial distributions, albeit with enhanced con- 
nectance. The incorporation of correlated inheritance has the effect of producing a longer 
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Figure 6: Plot of degree distribution for the 30000 unit resource system. Compared are the 
binomial distribution of the same connectance, C, and the binomial distribution expected for 
the null case with the space connectance C = 9 Q = 0.05 
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talled distribution that conforms closely to an exponential. 

A notable effect of the evolution is the dramatic increase in the mean degree compared 
to that of a random species set of the same diversity. There are a great deal more connec- 
tions between the extant species than is found in the null case. This is demonstrated in 
Fig. (6), where the evolved distribution can be compared to the binomial distribution that 
would be found for a random set of 29 species. The connectance values for all diversities 
are very much higher than general connectance of the interaction space 8q (see Fig. (8)) 
and they easily exceed the percolation thresholds that would be assigned to networks of 
these node/ diversity numbers (Albert & Barabasi, 2002). It is unclear whether greater 
stability is achieved through either lower than threshold connectance values or higher 
values (May, 1974) (McCann, 2000) (Haydon, 2000) (Dunne et al, 2002), but here the 
system naturally evolves towards the latter. It appears that the system has a strong ten- 
dency to evolve towards configurations that are highly mutualistic by assuming greater 
numbers of positive interactions. This is highlighted by the interaction strength distribu- 
tion that has clearly shifted from the null system normal with zero mean to a Gaussian 
form with a positive mean, Fig. (7). As we increase resource, and therefore the diversity, 
the mean connectance decreases as the system becomes progressively less likely to achieve 
large deviations from the null system expectation values, Fig. (8). In the limit of large 
diversity the connectance would of course conform to that of the interaction space, 9 . 



4 Discussion 

The concepts used in this model have been developed from original ideas presented in 
the Tangled Nature model (Christensen et al., 2002) (Hall et al., 2002) (Collobiano et al., 
2003) (Anderson & Jensen, 2005), where species interaction networks were constructed 
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Figure 7: Plot of the interaction strength distribution for the 30000 unit resource system. The 
evolved distribution deviates from the null normal towards positive strengths. 
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Figure 8: Plot of mean connectance, < 8 > against diversity. Each resource level is treated 
separately so multiple data points often occur for the same diversity value. Error bars represent 
the standard error. The lower dotted line marks the null system connectance, 9$ = 0.05 
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from dynamics instigated at the level of individuals as opposed to higher level group- 
ings such as species. Also, the use of a predefined, closed set of possible interactions 
between species allowed statistical comparisons to be made between evolved systems and 
null species sets. These principles have been incorporated in the new model along with 
modifications to improve the reality of the dynamics. 

The main modification is in the creation of the phenotype space. This has permit- 
ted a much larger set of distinct entity types to be available, 10 80 as opposed to 10 20 , 
and allowed us to introduce mathematically quantifiable correlations in the inheritance 
of interaction properties between parent and mutated offspring. Along with a more con- 
tinuous phenotype description, these attributes make the model more representative of a 
real biological system. 

The new model includes a generalised niche property in the form of an intra-specific 
competition term in the weight function. It was realised that this was necessary to allow a 
diversity to flourish without competitive exclusion reducing the system to a pair of highly 
mutualistic species. We make no attempt to describe any niche properties nor numbers 
of niches but the inclusion of this concept can be justified as a simplistic representation 
of a heterogeneous environment. 

With the above qualities in effect the primary aim of the research here is to investigate 
how the constraint of a conserved resource effects the dynamics and interaction properties 
of the species evolution. The mean species diversity exhibits a monotonic increase with 
respect to system resource availability which is expected for a real system with a hetero- 
geneous environment. This increase is also accompanied by a greater level of fluctuation 
which is consistent with a reduction in the degree of species permanence, as implicated by 
the lifetime distributions. Although the effect is small due to the inclusion of all mutant 
phenotypes in the distribution statistics, there is a clear decrease in mean lifetimes as 
we increase the resource level. To ensure this decline is not simply an artefact of the 
higher diversity systems producing more short-lived mutants, which are included in the 
statistics here, we repeated the analyses omitting contributions below various threshold 
lifetime values. In each case the decline was observed. We can infer from this that the 
decrease in mean lifetime is applicable to the more permanent members which is more 
representative of a real system. Species considered in the collection of field data would 
not be random mutants, but recognisable community members. If there exists a higher 
species turnover rate at higher resource levels then it is conceivable that the evolutionary 
rate is enhanced as well. There is evidence of a positive relationship between evolutionary 
rate and energy supply but the mechanisms suggested for its existence are not related 
to the concepts here (Pawar, 2005). It has recently been discovered though that higher 
diversity systems may incur a greater turnover rate of species and so a smaller duration 
probability (Emerson & Kolm, 2005). Our work certainly points to the same conclusion. 

The realised species interaction networks universally conform to exponential degree 
distributions which have also been reported in some real ecosystems. The long tails 
of these exponential distributions indicate the existence of small numbers of massively 
connected species which is highly unlikely to occur in the random binomial distribu- 
tions. Simulations performed with random mutations exhibit binomial degree distribu- 
tions, which implies that the correlations involved in inheritance are at the root of the 
exponential degree distributions. It is conceivable that the non-binomial distributions 
found in real ecological networks appear for similar reasons. The component species of 
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these networks, whilst not necessarily speciations of one and other, are likely to have 
traits or properties that could be considered as correlated. Invasive forms of community 
assembly would lead to systems with lower levels of correlation and so may exhibit less 
pronounced deviations from the binomial as a result. 

The high connectivity of the hub species contributes to the overall connectance of 
the network which evolves to statistically large values. It is unclear whether this level 
of connectance affects the stability and permanence of the networks as the system is 
continually attracted to achieving connections. It is interesting though that the mean 
connectance drops to a near constant level as higher resources/diversities are approached. 
This type of behaviour has been alluded to in field studies (Martinez, 1992) (Drossel & 
McKane, 2003) (Montoya & Sole, 2003) with a power law functional form suggested for 
the connectance-diversity relationship. A hyperbolic form has been suggested (McKane 
et al., 2000) for reasons of stability with deviations from the inverse power law being at- 
tributed to fluctuations via species invasions. Our model, due its inclusion of mutualistic 
interactions, gives an alternative explanation to this phenomenon. Our systems evolve in 
a manner that intrinsically construct networks with greater numbers of interactions. As 
diversities grow, the probability of achieving such high connectances diminishes leading 
to the observed relationship between the two properties. 

Making comparisons between food webs and the more encompassing species inter- 
action networks of our model is problematic though. The mutualistic interaction is an 
important feature here and is vital to the proposition that connectance is an evolution- 
arily favoured global property. Most field data is predator-prey based and so precludes 
those types of interaction. The studied networks may well be embedded within larger 
species interaction networks that involve mutualists and commensalists though so our 
reasoning still has validity. 



5 Appendix 

In order to create a predetermined and correlated interaction space the following process 
is used. We consider the interaction strength of a due to f3, J a/3 , to be determined 
additively through contributions associated with each of the L traits, 




The contributions J, will take values that are normally distributed so the prefactor 
-jt is used to ensure the final strengths are similarly distributed. 

To acquire these separate contributions we first couple each trait Tf to all others, 
from both a and (3 sets, to provide L modified trait variables, . This is achieved 
through the sum, 

L L 

T f = H b '« T ? + H c * T k rnod(WOOOO), (4) 

3=1 k=l 

where bij,Cik G {—1,1} are independent, randomly assigned with equal probability 
and held constant throughout. This coupling format allows any shift (mutation) in a sin- 
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gle trait variable to lead to a shift of the same magnitude in all modified trait variables. 
The coefficients, bij, Cik ensure that the resultant shifts in the modified trait variables are 
different for each initially mutated trait. The intention here is to map these modified 
trait variables to values that will represent the separate contributions to the interaction 
strength of Eq.(3). As will be demonstrated, the couplings assist in producing a quan- 
tifiable correlation, but they may also be justified as representing real interactions. The 
efficacy of a trait is dependent upon both its own, and the other interacting species, trait 
sets. For example, in pursuit of prey, the cheetah's high speed would facilitate its success 
but be dependent upon its stealthiness and the prey's capability of flight. The traits are 
interdependent and so lead to an interaction strength (or similarly, fitness) in a coupled 
manner. 

The separate contributions to the interaction strength of Eq.(3), are now determined 
by mapping each modified trait variable to a strength value, 

A: 7fU^,ie[l,16]. (5) 

For this we use L = 16 independent series indexed from to 99999 to allow a one to 
one mapping between the modified trait variable and the strength measure (see Fig. (9)). 
To incur correlations we use first-order auto-regression series (Hamilton, 1994) that are 
created using the stochastic iterative Markov process, 

x n+1 = px n + (6) 

where is a normally distributed random variable, p determines the correlation length 
of the series and nG[0, 99999]. This form of series has the overall property that its values 
conform to a normal distribution (once normalised) so providing trait-specific interaction 
strengths 3^ that are distributed in the same manner. 

The correlation function of the series has an exponential decay, 

C c = exp[-\n - m|/£], (7) 

where n, m are the indices of the series values. The correlation length, £ = 250 
is a metric of the index separation of the series values that is definable through the 
parameter, p from Eq.(6). With the use of full coupling of all traits to give a modified 
trait this function extends to separations in the phenotype space where the differences in 
trait coordinates sum to give an approximate two-point correlation function, 

~ exp [ Ll ^ i = 1 H (8) 

where, 

ATf = \b n (T? - Tf ) + b i2 (T? - 7?) + ... + b iL (T? - 7f)|, (9) 
and b^ are the parameters used in Eq.(4). 

If a mutation is imposed by the shift of a single trait variable then the numerator 
of Eq.(8) simply reduces to the length of that shift. Thus a mutated phenotype has a 
quantifiable correlation with the parent, with an exponential decay determined by the size 
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Figure 9: A section of a first order autoregression series. The values can be seen to have a 
degree of correlation. The series values are globally normally distributed so by use of the error 
function we can impose a threshold mark (straight line) which is a requirement for a coupling 
to be deemed to exist 



of the trait shift. This is where the trait couplings of Eq.(4) facilitate the quantification 
of the correlation. 

In a real ecological system, not all phenotypes interact directly with one another so 
only a certain proportion of the set of all interactions is allowed in the model. This is 
achieved by using the same process used to determine the interaction strengths except 
with some important modifications. Using a further set of autoregression series the pro- 
cess is repeated and the interaction is deemed to exist if, and only if, the two phenotypes 
have mapped values above a certain threshold, Fig. (9). As these values are normally 
distributed the error function can be used to determine this threshold so as to set the 
connectance of the interaction space 9q. For a particular phenotype, the proportion of all 
possible interactions it has across the phenotype space will be equal to this value, barring 
statistical fluctuation. The distribution of this set throughout the space is quite diverse 
with patches of varying size due to the random but correlated nature of the autoregression 
series and the interdependency of the traits. 

Phenotypic self-interaction is excluded in this model (intra-specific competition is con- 
sidered separately Eq.(l)) so the interaction set is manipulated to achieve this. When 
determining the existence of an interaction we use Eq. (4) as in the procedure used to de- 
termine the interaction strengths. But, there is a small but important modification. We 
can randomly set the parameter bij e {— 1, 1}, and set the second parameter = — b^ 
so that a pair of identical phenotypes will incur all modified trait coordinates to be zero. 
If we set our zeroth autoregression value low in comparison to the interaction threshold 
then self-interaction is entirely excluded. This effect extends beyond the phenotypic self- 
interaction to the species description as the correlations ensure that the exclusion occurs, 
with high probability at least, for phenotypes that are similar but distinct. The autore- 
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gression series are created with correlations running from the zero index point in both 
directions (ie. 0,1,2,... and 0,99999,99998,...), thus conforming to the periodic boundary 
of the trait variables. So, this species level negation of self-interaction is isotropic with 
regard to a mutation in either direction of a trait. 
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